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We estimate the energy density and the gluon distribution associated with the classical fields 
describing the early-time dynamics of the heavy-ion collisions. We first decompose the energy 
density into the momentum components exactly in the McLerran-Venugopalan model, with the use 
of the Wilson line correlators. Then we evolve the energy density with the free-field equation, which 
is justified by the dominance of the ultraviolet modes near the collision point. We also discuss the 
improvement with inclusion of nonlinear terms into the time evolution. Our numerical results at 
RHIC energy are fairly consistent with the empirical values. 



I. INTRODUCTION 

Color Glass Condensate (CGC) provides us with a 
theoretical foundation for the weak-coupHng descrip- 
tion of soft partons associated with highly energetic 
hadrons [J, The number of the soft partons grows 
with increasing energy due to quantum branching, which 
eventually leads to merging and saturation of the par- 
tons in the projectile at a certain energy characterized 
by the scale Qs- Such soft (wee) partons are relevant not 
only in the diffractive process in deep inelastic scatterings 
but also for thermalization in relativistic heavy-ion colli- 
sions . In the latter case, the dense transient system at 
mid-rapidity is initially created through the interactions 
between the soft partons from the incident nuclei. 

In the Au-Au collisions at the top energy y^s^Tw ^ 
200 GeV of Relativistic Heavy Ion Collider (RHIC), 
the relevant scale of Bjorken's x is roughly estimated 
as p_l/^s„„ ~ 10~^ with presuming that the bulk of 
the initial medium is composed of the gluons of mo- 
menta p_L < 1 GeV. Then, the phenomenological Golec- 
Biernat-Wiisthoff fit [1, Q i multiplied by a nuclear en- 
hancement factor A^/^ = 5.8 (see Refs. 0,(3) gives rise 
to — 2 GeV^. This implies that the bunch of the in- 
coming gluons with p± < Qs in the incident nuclei arc in 
the saturation regime. These gluons are to be described 
as the classical Weizsacker- Williams fields [1, Q in the 
first approximation. Such a classical field picture is ex- 
pected to be more reliable at the energy of Large Hadron 
Colhdcr (LHC), ^/s^ = 5500 GeV. At this energy the 
saturation scale for the Pb-Pb collisions will be around 
Ql ~ 5.2 GeV^ i.e., (1.6)^ times larger than that of 
RHIC. Since it will be shown that the energy density 
and multiplicity will scale with and Q^, respectively, 
the former will become 1.6'^ = 4.1 times larger and the 
latter 1.6^ = 2.6 times larger at LHC than at RHIC. 
At the same time, the typical time scale at LHC will be 
1.6^1 = 0.63 times shorter than at RHIC O, EH. In 
this paper we will address the initial stage of the nuclear 



collisions at RHIC energy taking the saturation scale as 
= 1-2 GeV^. [This uncertainty comes from the choice 
of relevant x.] 

There are a number of efforts to solve the classical 
Yang-Mills equations of motion in order to determine the 
early-time evolution in the heavy-ion collisions @, Q . The 
numerical studies have been performed quite successfully 
so far [13, [3, [IE, [3 and also the analytical tech- 
niques are developing recently [13, fTl- fl9t based on the 
McLcrran-Vcnugopalan (MV) model [2| in the CGC pic- 
ture. 

This paper is a continued attempt from Ref. [ist to 
give an analytical estimate for the energy deposit and the 
gluon production from the colliding coherent fields. It is 
well known that the initial gauge configuration contains 
only the longitudinal fields [1, H, [H, [13] ■ These longi- 
tudinal fields rapidly decay and the nonzero transverse 
fields are generated in the early-time evolution. How- 
ever, if we expand the evolution of the energy density in 
a power series of r from the collision time t = 0, we find 
that an ultraviolet (UV) divergence in each term, espe- 
cially in the zeroth order term, i.e. initial energy density. 
This divergence originates from the perturbative tail of 
the gluon distribution in the MV model. 

In Ref. [ll] an ansatz of the logarithmic form was pro- 
posed by one of the authors to resum this UV singular- 
ity. Here we take another way. Because the origin of 
the divergence is in the UV regime where couplings be- 
tween the classical fields are unimportant, one should be 
able to tame the divergent terms within the perturbative 
framework. This is indeed the case, as argued in Ref. [21 1 
in a different context and emphasized also in the numer- 
ical simulation [l6j . In this paper we shall resum all 
the most UV-singular terms, which results in solving the 
free-field equation. The solution is the Bessel function in 
the boost-invariant expanding geometry. At finite proper 
time T we can send the UV cutoff to infinity. Moreover 
we find that the energy density behaves as 1/r at large 
T, consistently with the free-streaming behavior. In fact. 
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the limit r reproduces the original UV divergence of 
the energy density as it should. 

Our present approach consists of the exact initial con- 
dition of the MV model and the following perturbative 
time evolution. The saturation effect is fully taken into 
account in the initial condition in this sense. We first 
neglect any nonlinear effect in the time dependence, and 
therefore the validity of this first treatment is limited to 
the very early time when the UV modes dominate the 
dynamics. This approach recovers the same analytical 
result as Ref. 2l| in a simple manner. Although the 
free-streaming region is no longer in the reliable range of 
the approximation, it is worth estimating the energy and 
the gluon multiplicity there and comparing them with 
the empirical values. Next we will examine the size of 
the nonlinear effect in the time evolution by including a 
mean-field type resummation in the equation. 

Our analyses should be useful to grasp deeper un- 
derstanding of the qualitative feature inherent to the 
Glasma [la], a transient stage from the coherent CGC 
state decaying to the thermalized plasma state. From 
our study two remarks are here in order. 1) Although 
the initial energy at the coUisionpoint in the MV model 
contains the UV divergence (l6l.[2]|. the momentum spec- 
trum of the energy content and the gluon distribution is 
well-defined as we will explicitly compute. The spectrum 
is quite informative on its own. 2) Another potential ap- 
plication is the analytical approach toward the Glasma 
instability p^ . [l^ . [2^ found in the numerical simula- 
tion [23| . Because the instability presumably resides in 
the very early time where our description of the time 
evolution works, it is a feasible strategy to examine the 
stability of field fluctuations on top of the perturbative 
time evolution. We will list more future outlooks in the 
final section. 

This paper is organized as follows: In the next section 
we introduce the MV model. Then, in Sec. IIIIl we cal- 
culate the correlation functions of the gluon fields and 
fix the MV model parameter corresponding to Qs for a 
given infrared regulator. Section HVl is devoted to a guide 
for our calculation procedures which are divided into the 
following four sections; the proper time expansion is dis- 
cussed in Sec. |Vl the energy density and the gluon dis- 
tribution at T = are evaluated in Sec. IVII the time 
evolution is convoluted in Sec. IVIIi and the time evolu- 
tion is augmented with nonlinear terms in Sec. lVIIll Our 
discussions and outlooks are in Sec. |IXl 



II. MODEL 

For the purpose of describing a longitudinally expand- 
ing system it is convenient to adopt the Bjorken coordi- 
nates of the proper time r and the space-time rapidity 77, 
defined respectively by 
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(1) 



It should be noted that 77 above is different from the 
pseudo-rapidity in momentum space which is often de- 
noted by r] in literature. The metric tensor associated 
with the Bjorken coordinates is Qtt = 1, Srp; = ^t'^, 
Qxx = Qyy = — 1, and zero otherwise. 

We will work in the radial gauge, At = 0, through- 
out this paper. The canonical momenta (chromo-electric 
fields) in this gauge read from the Lagrangian as [l^, [isl . 
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(3) 



Here we include r in front of C in Eqs. ([2]) and ([3]) origi- 
nating from \/\g\ in the integral measure. The following 
equations of motion are derived from Hamilton's equa- 
tions: 



drE' = - 
drE'^ = - 
with the Hamiltonian. 
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(4) 
(5) 

(6) 



These four equations ©-([S]) are the basic ingredients 
for the classical description valid in the early stage when 
small- a; partons are abundant and quantum corrections 
are still negligible; in momentum rapidity Y{= ln(l/a;)) 
space the validity region is bounded as \xi{\/as) <C F <C 
l/tts specifically. 

It has been well established that the initial condition 
is uniquely determined by boundary matching at the sin- 
gularities of the color source in the limit of vanishing lon- 
gitudinal extent due to the Lorentz contraction [1, 0, [2^ . 
The initial fields arc thus known as 
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where a^^^ and af^ are the pure gauge fields in the space- 
like region extending from the right-moving nucleus trav- 
eling on the axis and the left-moving nucleus traveling 
on the x~ axis, respectively. They are the gauge trans- 
formation from the light-cone solution [2^ : 
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(a; J.) = -—Voo{xA_)di V^{xa_) 



(8) 



with the Wilson lines, which are color SU(Afc) matrices 
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in the fundamental representation, defined by 



write down the expression, 



V^4x^)=V- exp 



dz 



— OQ 

+ 



(9) 



Here V± represents the path ordering with respect to 
x^. We remark that p'-^-' is a static color source in 
(i.e. x^-independent) and p'-^-' is static in x~ due 
to the Lorenz time dilatation. Although we assume 
p^^-'(a;^,a;~) cx 6{x~) and p^'^\x^,x'^) cx 6{x'^) in the 
end, we have to keep both the z-intcgral and the path 
ordering; which encompasses the random color distribu- 
tion along the longitudinal extent [J, [H [13, HI HI, HI • 

Now that we have fixed the initial condition for the 
equations of motion, we should have a unique solution. 
In principle the solution is given in terms of V and W 
through the initial condition, and any physical observ- 
able O at later time could be determined once p^^^ and 
p^^-* are given. Since the color sources p^^-* and p'-^-' will 
fluctuate randomly in event by event, observed O should 
be an ensemble average over the color source distribution 
>V[p(i),p(2)]^ that is. 
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This is a general form for the observable expectation 
value in the CGC framework. 

The MV model is founded on the Gaussian approxima- 
tion to W[p^^\ p'^'^''] as formulated in Rcf. [2^, where any 
averaged quantity is expressed in terms of the two-point 
function, 



(p(")"(a;^,z)p(")''(y^,z')) 

= .gV'(^) S'^^'Siz - z') S'^^Hx^-y^) 



(11) 



Here m and n distinguish the right-moving (1) and left- 
moving (2) nucleus, a and b are the color indices, and 
n'^iz) is the only dimensional parameter in the MV 
model. All the physical quantities are given in terms 
of (integrated) p^ and importantly p^ has a tight con- 
nection to the parton saturation scale Q^. 

With these preliminaries, we are ready to compute 
physical quantities of our interest in the MV model, once 
the model scale p^ and the coupling constant g are fixed. 
Before proceeding to the next section, let us elucidate a 
simple expression for a*^"^ for later use. Here we only 
deal with a^^^ of the right-moving nucleus and suppress 
the superscript in the rest of this section and in the next 
section. Similar formulae obviously hold for a'^' , too. 

Substituting Eq. Q into Eq. ([5]), we can explicitly 



ai ~ / dx Vj. 

^9 J-oo 



dip{x ) 



dx-'-^v.-nvi 



dx 



(12) 



F 1 



where Tp's are the SU(A*'c) algebra in the fundamental 
representation and Va denotes the Wilson line in the 
adjoint representation whose components are given by 
V^a"'' = 2tr[r^yr|yt] in terms of the fundamental Wil- 
son lines. 



III. GAUSSIAN AVERAGE 

In this section wc list the useful formulae for the Gaus- 
sian average over the color source distribution. Using 
Eq. ^2]) we can explicitly write the gauge field correlator 
in a form of 



dx-dy-^^'^"^'^^'^ y^^p"' '^y^^y ) 



dt 



X {vp.l^_\x^)vj.^;'\v^)) 

/oo 
dx- dtdy^L{x^-y^:)p^{x-) 
-OQ 

X Cadj(a;";a;^- yj^) , 
where we have defined 

' p°-{x^,x-) p^{y^,y-) 



and 



(13) 



(14) 



L{x^) = 



Xl = 



(27r)2 (fc^+TO2)2 



(15) 



The momentum integration in Eq. (jl5[) has infrared (IR) 
divergence. We regularize this by the gluon mass m ^ 
and regard m as a model parameter controlling the IR 
behavior. Then the integration can be performed analyt- 
ically. For convenience we introduce the notation [27l.l28t. 
T{x±) ~ 2L{0±) — 2L{xj_), which appears in the expres- 
sions for the physical quantities and is expressed as 



T{x^) = 



1 



Stt 



1 — m|a;j_| Ki{m\x±\) 
■\n[m^\x_L\^X^] , 



(16) 
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where Ki{x) is the modified Bcssel function, and A = 
g7E-i/2^2 ~ 0.54. The second fine is an approximate 
expression for small m|a;x|. 

Besides IR singularity, UV divergence also exists in 
d'^L{x±), {d^Y L{x \^ and so on, due to the Delta func- 
tion correlation 5^'^\x±^ in Eq. (fTT|) . as explicitly shown 
below in this section. 

In Eq. p^ . we have also defined the Wilson line cor- 
relator by 



(17) 



We note that Eq. is the building block and its it- 
erative use in the left-hand side of Eq. p?|) leads to 
C'adj(a^~; a;j_). The answer is already given in litera- 
ture [13, m, [11 as 



Cadj(a; ; = cxp 



(18) 

This factor encodes the effect of the multiple interac- 
tions with the color source in the nucleus. We see that 
Cadj{x^ ;x± 0±) ^ 1 as it should be due to color 
transparency. 

To simplify the discussion we assume here /i^(a;~) = 
d{x~)9{e — x~)^\/e with e ^ 0+ (this can be relaxed in 
more careful calculations) so that 



A- 



(19) 



dx fJ.'^ix ) 
Then the integration over x~ results in 

dx^ fj.'^{x^)Cadj{x^;x±) = /x^C'adj(a;±) (20) 

— oo 

with 



1 



exp 



am'^T{x±) 
We have used a dimensionless parameter 

_ A^c(gVA)^ 



am-Tix^) . (21) 



(22) 



in accord with Refs. [27|, |28|, |29|. Here in this paper, 
however, we regard the parameter a as the IR parameter 
controlling how much of the non-perturbative effect in 
the initial fields contributes to the momentum integration 
once the saturation scale is fixed [see discussions below] . 
Finally we find from Eq. ([T^ . 



anx±)a'^{y^) 



5'''g^^i\ C^i{x^-y^) d^d^Lix^- y, 



(23) 



Let us discuss a rough estimate for the parameter a 
here. We will check at the end of this section that g^iJ,A is 
nearly identified with the saturation scale Qs (which also 




FIG. 1: Cadj(a^x) as a function of dimensionless m^/a\x±\ = 
^Nc/2g'^^iA\x±\ in the case that a = N^ig'^ ^ia^ / {2m^) is 25 
which is the same choice as Ref. [2^. The solid curve is the 
result from the numerical integration and the dashed curve 
represents the approximation by Eq. (|16|l which eventually 
blows up in the IR region. 



depends on a). The IR cutoff is difficult to fix uniquely 
but a natural expectation is to ~ g^j^A^ from the compar- 
ison between our cutoff scheme and the color neutrality 
as a result of the quantum evolution p6l . [sij with an un- 
certain prefactor. The strong coupling constant at the 
scale Qs at RHIC energy is as — 0.3 (or g ~ 2). This 
implies a ~ 3 • 5^/2 ~ 6. The physical meaning of this m 
value is that the classical description breaks down there. 
On the other hand, if we take to of order of the con- 
finement scale ~ 1 fm, i.e., m ~ 200 MeV, we would 
have a ~ (3 • 1 - 2 GeV^)l{2 ■ 0.2^ GeV^) = 37-75 
for = 1 ~ 2 GeV^. In this estimate m is the scale 
where the CGC picture in the perturbative regime breaks 
down. From these considerations the reasonable range 
for a would be 10-100 at the RHIC energy. It should be 
noted that a used in the numerical simulation may be 
much larger; there to is provided by the system size of 
order 10 fm and a may be then hundreds times larger. 
Nevertheless, because the IR property in the numerical 
implementation of the MV model is totally different from 
the analytical formulation [2^ . we cannot make a direct 
comparison to the numerical results. In this paper we 
will specifically choose four cases: a = 10, 25, 100, and 
500 for comparison, keeping in mind that the physical 
value is around a = 10 ~ 40. 

It is important to note that the above Cadj {x\_) exactly 
corresponds to the gluon propagator [2^ and should be 
distinguished from 



Ca{x±) = exp[— aTO^r(a;j_)] 



(24) 



which appears in the gluon distribution [27|, |28|, [S^] . In 
the case of dense-dilute collisions it is possible to convert 
Ca.di{x±) into Ca{x±) by gauge rotation but this is not 
the case if the dense-dense collisions are concerned in a 
symmetric way. 
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The part involving C'adj(a'j_) describes the transverse 
color correlation in the random distribution. It shows 
moderate damping in a distant region as plotted by the 
(red) solid curve in Fig. [T]for a = 25. We show the ap- 
proximation using Eq. (|16p by the (green) dashed curve, 
and Ca{x±) by the (blue) dotted curve for reference. 
From the comparison between the solid and dotted curves 
we immediately realize that C'adj(a;±) retains a long- 
range tail than Ca{x±). Indeed, we can readily notice 
that r(a;^) — > {2'Kvn?)~^ for m|a;^| » 1 and thus 



Cadj(a;_L) 



27r 



1 - e--^/^-) 



(25) 



which is nearly 0.50, 0.25, 0.063, and 0.013 for a = 10, 
25, 100, and 500, respectively; the asymptotic value is 
fixed by the IR cutoff a. 

Equation p3p is a useful expression. Here we enumer- 
ate quantities necessary for the evaluation of the initial 
energy; the equal-point gauge field correlator is 



1 



and the derivative correlations read 
1 



(26) 



N,g^fi\6,,6kiid^LiO^)f 



(27) 



and also 



anx^)^'a%x^))^--S-'g'^,'^ 



N,g^^i\5^,{^''L{Q^)f 



25,j{dyL{Oi_) 



(28) 



where we can explicitly estimate the integrals as 
<Pk I 1 



(2^)2 fc2 



:^ln-, (29) 
27r m 



and 



1 



= 

(27r)2 47r 



(30) 



for A ^ m. Wc see that the derivatives of the correlation 
L{x±) at the origin contain the UV divergence. 

Finally in this section, let us make the connection be- 
tween the scale parameter, g^fiA, and the saturation scale 
Qs- The definition of the saturation scale is not unique. 
We apply the same scheme as adopted in Ref. dH us- 
ing the Fourier transform of Ca{x±), which we denote 
as CA{k±). Since fc^Cyi(fc±) is interpreted as a gluon 
distribution function, we define the saturation scale as 
the peak position of this function. Bear in mind that we 




FIG. 2: k^CA(k±) as a function of the dimensionless variable 
k± / (m^/a) where my/a = \/ Nc/2 g^fiA for a — 10, 25, 100, 
and 500. 



a 


m 




10 


0.64Q, 


1.650, 


25 


0.36Q, 


1.46Q, 


100 


0.14Q, 


1.13Q, 


500 


0.050Qs 


0.90Q, 



TABLE I: IR cutoff m and the MV model parameter 
determined from k\CA{ki_) for various a. 



have to keep the same definition for Qs in order to make 
a consistent comparison with other empirical analyses. 
[Qs itself is not gauge invariant.] 

Figure [2] is a plot for k\CA{kA_) as a function of 
k±/{my/a) where my/a = ^JN/J2 g^ \i a- The peak po- 
sition defines the saturation scale. From Fig. [2] it reads 
g,,/(mV^) = 0.497, 0.558, 0.720, and 0.903 for a = 10, 
25, 100, and 500, which correspond to m = 0.64(5s, 
0.36(3s, 0.14(3s, and 0.050(3s, respectively. We summa- 
rize these relations in Table HI 



IV. CALCULATION STEPS 

In the subsequent sections we will elucidate the time 
evolution of the initial energy density and the gluon dis- 
tribution. Since the computation process is involved, we 
outline the derivation in advance here. 

Step I) We will carry out the proper time expansion; 
the equations of motion arc solved in a power series of r. 

Step II) We will take the Gaussian average over the 
color source distribution order by order of t. We will 
then find that each term contains the UV singularity. 
We decompose the divergent initial energy density into 
the Fourier modes which give the gluon distribution if 
divided by the gluon energy. We take full account of the 
saturation effect. 

Step III) Wc will figure out the time evolution for the 
UV modes from the equations of motion. We confirm 
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that the UV singularity is regularized at finite r. and fix 
the initial condition for the time evolution by matching 
it to Step II). 

Step IV) We will improve the time evolution by in- 
clusion of the nonlinear terms in the Gaussian approx- 
imation. We finally read the initial energy density and 
the gluon number at the formation time when the time 
dependence shows free-streaming behavior. 



Then we can explicitly write down the second-order 
contributions to the initial chromo-electric and chromo- 
magnetic fields: 



and 



(2) 



-^(2) - 5^J(0)-P'ji)(2) 



i^j(0)^j(0)-E^(0) ' 



(36) 



V. PROPER TIME EXPANSION - STEP I) 



B\2) 
^(2) 



jri{2) - e'H^j(0)£'(o) : 



12(2) 



3-Oj(0)^i(0)S('o) 



(37) 



We shall begin with a naive expansion for the classical 
fields in terms of the proper time r as in Ref. [l3], and 
in the next section we resum the leading order terms in 
the UV singularity coming from the Gaussian average. It 
should be stressed that the gauge fields for a given color 
source are free from the singularity and have no difficulty 
in the r-expansion. 

The first order terms are the contributions at t = 0, 
that is, the initial fields as given in Eq. ([7]). Let us 
now evaluate the associated chromo-electric and chromo- 
magnetic fields. It is trivial to read the initial chromo- 
electric fields from the initial condition as 



^(0) - ' 



^(0) ~ 



and noting a^"'''s are pure gauge solutions, we find the 
chromo-magnetic fields being 



J2)- 



(31) 



-B(o) - , Bj^-^ = i^i2(o) 



(32) 



where 



,0)=-..g([aP^af] + ["f^«^])• (33) 



These results are recognized as the fundamental prop- 
erty of the Glasma initial state; the longitudinal fields 
between two color sheets are predominant in the initial 
state of matter. 

To proceed further away from r = 0, we apply the 
expansion, i.e.. 



0{t) = Y.O 

n=0 



in)'' 



(34) 



for arbitrary fields O given in terms of the classical gauge 
fields. It is easy to confirm that the terms for odd n are all 
vanishing due to time reversal symmetry. The non-trivial 
contribution to the gauge field starts at the second-order 
terms which turn out to be 



■ 

A(2) 



1 Tpi 
2-^(2) ■ 

2^(0) 



(35) 



Here we have defined the anti-symmetric tensor e^^ ~ 
— e^^ = 1 in the transverse coordinates. The duality rela- 
tion between the electric and magnetic fields is manifest 
itself in Eqs. (jM)) and ([57]) . 



VI. INITIAL ENERGY DENSITY AND GLUON 
DISTRIBUTION - STEP II) 

The energy density is an ensemble average of the 
Hamiltonian ©, i.e. e = {TL). In this section we compute 
e and decompose it into the Fourier components. The 
merit of the Fourier decomposition is that each Fourier 
mode is free from UV singularity and we realize how the 
fc^-integration of the spectrum diverges. Besides, we can 
define the gluon distribution associated with the field in- 
tensity at each momentum. 



A. Zeroth Order 

Right after the collision at t = only the longitudi- 
nal fields along the direction have nonzero values, and 
therefore we calculate the following energy density: 



^(0) 



^(0)^(0)J 



(38) 



Noting that r(a;j_— y^) ^ in the limit of x^, 
we can drop C'adj(a;^ — yj^) in Eq. ((38|) . With a short- 



hand notation as (a,-"'''°a("''^^ 
express the energy density in a form of 



) = ,5'""(5°'',5y (aa), we can 



MO) 



2iVc(iV2 



5V1 



In 



A 



, (39) 



from Eq. ^ for = 3. 

It is interesting and practically useful to look into the 
energy content in momentum space, which is also nec- 
essary when we consider the gluon distribution. To this 
end, let us calculate the following quantity: 



£(o)(fcj-) = :^(tr[i?['„)(-fcx)£;("o)(fc±) 



(40) 
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so that the fcj^-integration of e(o)(fc±) recovers the en- 
ergy density £(o)- Then we have to evaluate the spatiaUy 
separated correlation function of the chromo-clcctric and 
chromo-magnetic fields: 



trE^^,^{x^)El,iy^)) = -N^iN^ - l)g 



"(0) 



\dlL{x^ - y^)f + 2{did2L{x^ -y^)f 



(41) 



+ {dlL{x^-y^)Y 

Computing the chromo-magnetic part in the same way, 
we obtain almost the same expression as Eq. (PT|) but 
with the opposite sign in front of the [did2L{x± — y ^y\^ 
term. This term cancels out in the sum of the chromo- 
electric and chromo-magnetic parts. The Fourier trans- 
formation of the sum gives rise to 

<Pqil_ d^Q2± d^Q3± C'adj(gi±)C'adj(g2±) 



(2^)2 (2^)2 (2^)2 (q2^+m2)[(fcl-q3^)2+: 



(42) 



We denote k'^ ^ k± Qi±— Q2±^ ^^'^ ^adj(g±) is the 
Fourier transform of Cadj(3;j_). Wc can further perform 
the q3j^-integration to reach the following: 



£(0)(fci) = ^^c(iVc'-l)5Vi 



(fq^^ (fq2± 
(27r)2 (27r)2 



^ Cadj(gix)C'adj(g2±) 



+ 4m2 -I- k'. 



^k'l + 4m2 - k\ 



(43) 



At this point it is easy to confirm that Eq. (|43|) repro- 
duces Eq. ((39|) when integrated over k±. In the presence 
of the fc^-intcgration. we can change the variable from 
fcj^ to Aj'l = fc^ — <7]^^ — and then wc can perform the 
1i 2±-iiitegrations independently of the fcx-integration. 
Using the sum rule, 



(27r)2 



■ adj(gi,2±) = adj(a;± = Oi) = 1 , (44) 



we have 



1 

47r 
1 

8^ 



7Ve(^c'-l)5VA 



In- 



(Pkx_ 1 
(27r)2 

A'' 



pcrt{k±/m) 
(45) 




k\ I m a 



FIG. 3: mY^fcxCadj(fe±) as a function of the dimensionless 
variable k±/{m^/a) for a — 10, 25, 100, and 500. 



where we have used A 3> m and this is exactly identi- 
cal to Eq. ([55]) . Here, we have defined the perturbative 
contribution by a function. 



1 



:ln 



(46) 



It is worth mentioning that T^crtlC ^ 0) ^ 0.5 which is 
finite owing to IR rcgularization. 

Wc need to evaluate Cadj(fc±) appearing in the inte- 
gral P5|) . In Fig. [3] we plot m^/ak±Ca,dj{k±) as a func- 
tion of kj_/{m^/a). As is already mentioned previously, 
Ciidj{x±) has a constant tail (27r/a)(l — exp[— a/(27r)]) 
at large distances due to the IR cutoff, which gives rise 
to a Dirac delta function S^'^^{kj_) in momentum space. 
C'adj(fcj_) represents the effect of the multiple scattering. 
For smaller a with fixed g^i^A, larger part of IR region is 
cutoff and therefore Cadj {k±) spreads over smaller k± re- 
gion and has larger weight of 5^2) (fcj_). For larger a, more 
multiple interactions are effective and C'adj(fci) has the 
larger spread and smaller weig ht of (5(2)(fcj_). The per- 
turbative limit corresponds to no multiple scattering. 

Using the numerical results shown in Fig. [3] and 
the Dirac delta function with the weight (27r/a)(l — 
exp[— a/(27r)]), we can explicitly evaluate e(o)(^-L) Per- 
forming the numerical integration of g^j. ^^^^ Q2±- the 
perturbative limit where C'adj(gi2j.) factor is replaced 
with (27r)2(5(2)(q^^^) in Eq. ([33]), the fcj_-spectrum of the 
energy content is simply fixed by Tpe,-t{k±/m). Wc show 
this perturbative limit of Tpci-tik^/m) by the (black) dot- 
ted curve in Fig. [H 

The saturation effect changes the functional form of 
'^ort('f) into the one depending on a which we denote by 
T(a;^), i.e.. 



£{o)(^±) 



1 



47rm2 



N,{N^-l)g^fi\Tia;k^/m). (47) 



Figured shows T (10; by the (red) solid curve, T (25; ^) 
by the (green) dashed curve, T(100;^) by the (blue) 
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order in the fields. That is, 



10 12 14 16 18 



FIG. 4: T(a; ^) in the initial energy density as a function 
of the dimensionless variable k±/m); the (black) dot- 
ted curve represents the perturbative limit TpcvtiO (where 
C'adj('7i 2±) is replaced with {2tv)'^ S^^\q-^ 2±))' ^^^^ (red) solid 
curve for a — 10, the (green) dashed curve for a = 25, the 
(blue) dotted curve for a = 100, and the (purple) dot-dashed 
curve for a — 500. 



dotted curve, and T(500; ^) by the (purple) dot-dashed 
curve. We notice at a glance at Fig. [4] that the saturation 
effect is larger for larger a in the low kj_ region. Although 
the functional form of Cadj(fe±) has a large suppression 
for small a in view of Fig. [3l the resulting spectrum in 
Fig.[3]turns out to be closer to the perturbative one when 
a is small. This is because Cadjl^i 2_l) contains the Dirac 
delta function S^'^^q^ 2±) with a larger weight for smaller 
a. 

If we integrate the above e(o)(fc_L) over k± the initial 
energy density becomes a-independent and should be re- 
duced to Eq. (|45l) . As we will discuss later, the spectrum 
should follow the perturbative time evolution, and then 
we will have the energy density as a function of time that 
is a-dependent and UV finite (but m dependent). 



B. Second Order 

We are now going into the contributions of ©(r^). 
There arise two kinds of terms with and without spatial 
derivatives. The derivative terms bring about the UV 
singularity from the Gaussian average in the MV model. 
The non-derivative terms have more gauge fields in place 
of derivatives. We will discuss each contribution in order. 



1. Longitudinal fields 

We calculate first the longitudinal field energy of 
0{t^). Because the longitudinal field has an 0{t^) con- 
tribution, the second-order terms in the energy come 
from the cross terms between the zeroth and the second 



pi 



tr[2i?,"o)i?r2)+25(o)5?2)] 



'(0). 



(48) 



The derivative part yields 



= 2g'^N,{N^ - l){aa) (ad^a) 



.6,4 



In- 



A 



(49) 



The non-derivative part, on the other hand, gives 
1 



^(0) K(0)' [^j(0).£^(''o)]] 

-^(o)K(o)J^.(o),i?;'o)]] 



InA 

m 



(50) 



The sum of these gives rise to the second-order term of 
the longitudinal field energy: 



27r^ TO Stt'' 



In- 



A 



(51) 



2. Transverse fields 

The computation for the transverse fields is parallel to 
the previous case. Since the Hamiltonian ^ possesses 
l/r^ in front of the transverse fields, the square of the 
0{t'^) fields contributes to the energy at O(t^). Note 
that there is no 0(t") transverse field. We thus evaluate. 



'-(■2) 



(2)^(2) 



^(2)^(2)] 



-Oj(0)^(0)^i(0)^[o)] 



(52) 



The derivative terms give 

\itr[d,El^d,El^+d,Bl^d,Bl^] 
= g^N^{Nl - l){aa){djadja) 



5V^T^A^hrA_gio^6^. 9 

47r^ TO OTT'^ 



In- 



A 



(53) 



which is nothing but a negative half of the corresponding 
derivative part in the longitudinal fields. We can also 
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check that the non-derivative terms are a negative half 
of the longitudinal counterparts as well. Hence we find 



-'{2) 



2^(2) 



(54) 



VII. TIME EVOLUTION - STEP III) 

We must take the limit of A ^ cx) in the end because 
the UV singularity is not physical in contrast to the IR 
property which is physical. In the situation with large 
A only the derivative terms oc h? are predominant in 
the O(t^) corrections. We shall therefore drop the non- 
derivative contributions for the moment. In fact we could 
take account of those effects and find them negligible in 
the early-time region. 



A. Energy Density 

In Ref. [l^ one of the authors proposed the following 
logarithmic ansatz to resum the UV diverging terms: 



e(r) 



^(0) 



-(2)T 

3 

47r2 m 



A2 

2 



In ^ I In ^ - -AV^ 



A_^ 1 

A2 



i2 + im^A^T^ 



A nice feature of this ansatz is that the A 
well defined, i.e., 



1/2 n4 3 
9 



In- 



• (55) 
oo limit is 

(56) 



This is a concise pocket formula. We will later confirm 
that Eq. (|56p works well to give a nice approximation to 
the perturbative time evolution of the Bessel function. 

In the rest of this section let us elaborate the correct 
resummation of the the highest-order divergent terms at 
each order of r. The idea is simple: We have under- 
stood that the UV singularity is attributed to the spatial 
derivative. Then, it is natural to anticipate that the so- 
lution of the equations of motion with the most singular 
spatial derivatives retained resums the UV singularity. 
This is indeed the case, as is explicitly done here. 

The highest-derivative contribution is the solution of 
the equations of motion: 



T 



(57) 
(58) 



From the above differential equations for Ai and A^j , the 
time dependence is deduced as 



Ai{T,kj_) = Ai(o) Jo(fc_LT) , 
A^{T,k±) = A,j(^2)-^Ji{k±T) . 



(59) 



We have dropped the transverse projection because 
the initial condition £'^'2) given as Eq. p6|) is transverse 
in the leading order of the spatial derivative. We can find 
E^ and E^ with the use of the definition of the canonical 
momenta ([2|) and Q. In the leading order of the spatial 
derivative we have = diA2 — 92^1 and = e^WjA^-j, 
so that B"^ and B^ have the time-dependence same as Ai 
and Ari, respectively. In summary, we have 

E\T,ki_)^ElAk^)Uk^T), 



2t 



E'iT,k^) ^ E'{k^)^Jl{k^T) , 

B'^{T,ki_) = B'i^^{ki:)Mk^T) , 

B\T,k^) ^ Bl2.{k^)^Mk^T) . 



(60) 



It follows that the longitudinal field energy component 
can be summed up to be 

£^(r,fc^) - -^g^^\Tia;k^/m)[Jo{k^r)]^ (61) 

with A^c = 3 substituted. It is easy to see that the 0{t^) 
term exactly corresponds to the A^ terms in e^2) if 
panded. Likewise, the transverse field energy becomes 



6 



,gV^r(a;fci/TO)[Ji(fcir)]', (62) 



\r,k^) 



whose coefficient is fixed so as to reproduce Eq. ((54| in 
the T-expansion. The total energy density at each mo- 
mentum thus evolves as 

e(T,fci) = e^o){k±)[[Mk±T)]^ + [Ji(fc±T)]') , (63) 

and the integrated one is 



e(r) = -!- [ dkj_kj_e{T,k±) 
Ztt J 

= ■ -^{g'^t^A)'^ lE{a;mT) , 
where we have defined a dimensionless function. 



(64) 



POO 

lE{a;mT)= / d(^T{a;0{[M^mT)]'+[j,{^mT)]'} 



(65) 

Although our derivation is rather straightforward, this 
time dependence according to the Bessel function is ex- 
actly the same as the one given in Ref. [2l|. It should 
be noted here that this Bessel function form has been 
discussed in the earliest work in R ef. [811 and confirmed 
in the numerical simulation in Ref. |12l | in the weak field 
limit. However, we do not assume the weak field limit at 
all in our argument but the UV dominance at r ^ l/Qs 
is sufficient to justify this time dependence. The same 
evolution is also found in the Abelian model in Ref. [l^ . 

Our expression for the initial energy density (j64p is a 
function of the proper time r, the QCD coupling g, the 
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FIG. 5: [Jo(a;)]^ + [Ji (a::)]^ as a function of a:: which goes to zero 
as 2/(tix) at large x. The dotted curve starting from 1 rep- 
resents [Jo(3;)]^ and the other one starting from represents 



MV model scale g'^fiA, and the IR cutoff m or its ratio a 
to g'^fXA- The Bessel function [Joik±T)\^ + [Ji(fc_LT)]^ en- 
compasses the time dependence which we plot in Fig. [5l 
Interestingly enough, the sum of two is a smooth func- 
tion, though each of [Ja{k±T)]'^ and [Ji{k±T)]'^ is an oscil- 
latory function of kj_T. This is reminiscent of the trigono- 
metric function. Moreover, we can prove that 



jk I 



(66) 



asymptotically in the fcj^r oo limit. This analytical 
feature results in the behavior of e(T — > oo) oc 1/t that 
is the same as the free-streaming expansion [2l| . 



B. Gluon Distribution 

As for the gluon distribution, we shall adopt here a 
working definition in accord with the preceding works 0, 
H, [13, Elj Hi nil- That is, we assume a harmonic os- 
cillator and count the number of quanta. The concrete 
procedure is that we divide the field energy with momen- 
tum fcj^ by the (free) gluon energy quanta k±_ to infer the 
gluon distribution contained in the fields, which is consis- 
tent with the multiplicity computation by the reduction 
formula [13] ■ We define in this way the gluon distribution 
as 



n(r, k±_) 



1 



:e(r, fc_L) , 



(67) 



where we put the gluon mass m. From this we readily 
reach the expression for the gluon distribution as 



n(T) 



1 

3 



dk±k± n{T, kj_) 



(68) 



FIG. 6; {mT)lE{mT) and {m,T)lN{mT) used in Eqs. (|70|l and 
(|7ip . respectively, as a function of m^/ar for various a. 



where 



/Ar(a; mr) 



(69) 



It should be mentioned that dn/cfik at r = goes like 
cx (1/fc^) ln(fc_L/m) for large k± which reproduces the 
perturbative results known as the Gunion-Bcrtsch multi- 
plicity [^, [35| . To compare with the observed multiplic- 
ity, we have to assume the parton-hadron duality and the 
entropy conservation to interpret the initial gluon distri- 
bution as converted to the particle multiplicity. 



C. Numerics 

We are now ready to give a numerical estimate for the 
initial energy and the multiplicity from our formulation. 
The key equations are Eqs. dMj), dM]), and 
From these expressions we can write the total energy and 
the multiplicity per rapidity as 



dE{T) 

dr/ 



dN{T) 

drj 



:{g'^HA)^{mT)lEia;mT) , (70) 



ttR\ Tn{T) 



^[g'^^iA)^{mT)lN{a]mT) . (71) 
TT^m^ g'^ 

We plot (mT) I E (a; mr) and {mT)lN{a;mT) appearing 
in the above expressions in Fig. [6] for various a us- 
ing T{a; k±/m) shown in Fig. U) It is clear that both 
(rnT)lE(mT) and {mT)Ii^{mT) behave approaching con- 
stant for large mr as we have mentioned before. 

The validity of the CGC (classical field) approxima- 
tion breaks down when the system becomes dilute near 
g^fiAT ^ my/ar ^ 1. In the very early stage of the evo- 
lution, on the other hand, the UV modes dominate in 
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100 



10 



0.1 



0.2 



0.4 



0.6 



a 


dE/d-q (GeV) 


dN/dri 


c 


10 


(4.9 - 14) X 10^ 


(2.5-5.0) X 10^ 


3.04 


25 


(4.9 - 14) X 10^ 


(4.0 - 7.9) X 10^ 


4.86 


100 


(3.5 - 10) X 10=* 


(5.6 - 11) X 10^ 


6.88 


500 


(2.6 - 7.3) X 10^ 


(7.7 - 16) X 10^ 


9.47 



TABLE II: Total energy and gluon distribution estimated in 
our formulation for various a with the choice of — 1- 
2 GeV^. [Lower and upper values correspond to = 1 GeV'^ 
and Ol = 2 GeV^, respectively.] We also list the gluon liber- 
ation factor c. 



VIII. IMPROVEMENT - STEP IV) 



FIG. 7; Comparison between the log-ansatz (ln[2/mr])^ 
(green dashed curve) and the numerically evaluated J_E(mr) 
(red solid curves) for a = 10, 25, 100, and 500 from the top 
to the bottom. 



the physical quantities and the logarithmic ansatz (|55[) 
viforks quite well. We see this fact from the comparison in 
Fig. [7] which is the same plot as Fig. [H] without the factor 
mr. This is quite a non-trivial statement. We started 
from r = with the exact initial condition provided by 
the MY model. Then we found that, as long as ra^/ar 
is much smaller than unity (or simply mr < 1), only the 
UV modes control the time evolution of the energy and 
there is only little effect from the IR sector which changes 
with a. Hence, the remaining a dependence is only the 
trivial scaling in terms of m with ra^fa = Nc/lg^ ^jla 
fixed as long as r is small. 

From Fig. [6] we can give physical numbers for the en- 
ergy and the multiplicity deduced in the l/r-scaling re- 
gion. We first set itR\ « 150 fm^ and g ~ 2 which are 
less ambiguous. 

In the case of a = 10, we have obtained m ~ 0.64Qs 
and (g^fiA)^ — ^.lOQ^. The asymptotic values of 
(mT)lEio,,mT) and {mT)I]\r(a,mT) at t oo are re- 
spectively 1.43 and 0.47. Thus, we have dE/drj = 
1.4 X 10^ GeV and dN/di] = 5.0 x 10^ in the case of 
Q2 = 2 GeV^ IfQl = 1 GeV^ dE/drj and dN/drj get 
smaller by factors 2^^^ and 2, respectively. For other 
values of a, we can get the estimates in the same way. 
We summarize our results in Table |TT1 We note that c 
listed in Table |TT] is the gluon liberation factor defined 
by MM 



1 dN 



drj 



(72) 



It is interesting to note that the energy is less sensitive 
to a in the range of our interest while the multiplicity 
rises significantly with increasing a. This is because a 
controls the IR scale m to which the average gluon energy 
(dE / drj) / (dN / dr]) is proportional. In other words, with 
increasing a we include more IR modes, which contribute 
to dN/di] more than dE/drj. 



Let us consider here a possible improvement of our es- 
timate. The UV dominance cannot last at later time, and 
the nonlinear effect may become manifest in the interme- 
diate stage of the evolution. We then need include less 
singular terms from the nonlinear interactions appearing 
in Eq. ([?T|) . In order to estimate the size of this nonlin- 
ear effect we take the following strategy: In the Gaus- 
sian approximation we adopt here, the nonlinear terms 
in the Yang-Mills equations of motion produce the mean- 
field contributions like the mass term, that would modify 
Eqs. dSZl) and ^ as 



(73) 



drE^' = t(02 - K^g^N,{aa))Af , 



d,E^ = -(^' - K^g^N,{aa))A^ , (74) 



where re^ and the coefficients coming from con- 

tractions of the gauge fields; AfAjA';^ S"'^Sij{aa)A^ 
and A^A^'jAj" (5°^%(aa)A^'=, etc. Noting that the 
initial Ai has no correlation with Af in the leading order 
of the derivative, on the one hand, we find = = 2 
by these replacements in the equations of motion, though 
this might underestimate the effect. On the other hand, 
if we determine and to reproduce the expanded 
energy in Eqs. ([5T|) and ([5i)) up to 0{t^), we conclude 
= 7, which we shall adopt here. 

In Eqs. ([75)1 and ([7^ . we have estimated (AiAi) at 
T = and we did not include (A^^,,) which is 0{t^). 
From Eqs. and (gSl) we know that g'^Ndaa) ^ 

(7Vc/47r)(.g2/i^)2ln[A/m] = (TO2a)/(27r) ln[A/m] with a 
introduced in Eq. ([22|) . Our discussion on the log-ansatz 
suggests that this ln[A/TO] is taken over by ln[2/mT] at fi- 
nite r (see Fig. [7]), and it takes the value around In 20 ~ 3 
to In 200 ~ 5.3 for mr = 0.01 ^ 0.1 of our interest (so 
that my/ar is below unity). During this interval, there- 
fore, we may replace this logarithmic function by a con- 
stant of order 0(1) ~ 0(10), in effect. Let us assume 
rather arbitrarily the value tt for ln[A/m], which is within 
the range from 3 to 5.3, so that the mean- field effect ap- 
pears only through fc^ — *■ fc^ -|- 3.5m^a in the argument 
of the Bessel function. 

Figure [5] shows how (TOr)/£;(a; mr) changes due to the 
replacement of /c^ — s- /c^ + 3.5m^a. We summarize our 
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FIG. 8; Comparison between the log-ansatz (mr){ln[2/mr])'^ 
(green dashed curve) and the numerically evaluated 
{mT)lE{a', mr) with (red solid curves) and without (black dot- 
ted curves) k'i ^ kl + S.Sm^a for a = 10, 25, 100, and 500 
from the top to the bottom. 



a 


dE/dr] (GeV) 


dN/dri 


c 


10 


(2.8 - 7.8) X 10^ 


(0.7-1.4) X 10^ 


0.86 


25 


(2.5 - 7.0) X 10^ 


(0.9 - 1.8) X 10=* 


1.10 


100 


(1.6-4.5) X 10^ 


(1.0 - 1.9) X 10^ 


1.18 


500 


(1.0 - 2.9) X 10^ 


(1.0 - 1.9) X 10^ 


1.17 



TABLE III: Total energy and gluon distribution estimated in 
our formulation with nonlinear terms for various a with the 
choice of Q'i = 1-2 GeV^ 



final results with nonlinear terms in Tab. lIIIl It turns out 
that the nonlinear effect leads to an interesting feature 
that the dependence of dN/drj on changing a is mild, 
while dE/drj retains almost the same (or stronger) de- 
creasing behavior as it is previously in Tab. |TT1 



IX. DISCUSSIONS AND OUTLOOKS 

We have shown by explicit calculations in the MV 
model that the UV components are dominant in the en- 
ergy density at t = produced by the collisions. From 
this observation we have evolved the exact initial energy 
provided by the MV model via the free field equation. 
We did not require the assumption of small field ampli- 
tude, but the UV dominance allows us to drop nonlinear 
terms in the early time evolution. 

To extend the validity region of the estimate, we have 
augmented the time evolution with the effect of nonlinear 
terms in the Gaussian-type approximation. Assuming 
that the IR cutoff m arises from the quantum evolution 
m ~ g^A or from non-perturbative dynamics m ~ Aqcd , 
we have the IR parameter a = 10-40. In this range of a, 
we have found that the nonlinear terms tend to reduce 
the a-dependence of the physical observables. Let us then 
take a = 25. Regarding the relevant saturation scale. 



we take = 1 GeV^ here because the discussion in 
Ref. [m suggests that Ql is close to 1 GeV^ than 2 GeV^ 
Thus, what we get for a realistic situation at RHIC is 
dE/drj ~ 2.5 X 10^ GeV and dN/drj ~ 0.9 x 10^ with 
c = 1.10. Our result favors a larger value for c than the 
early numerical simulations [l3|, but is consistent with 
the analytical calculation c w 21n2 w 1.4 37 1 and agrees 
well with the latest numerical evaluation 11 , [s^l ■ 

Our calculation seems to underestimate the total mul- 
tiplicity whose empirical value is ~ 1150 at RHIC. We 
should, however, be aware that dN/drj here signifies the 
total number of gluons produced initially. Therefore, 
though dN/drj should not be far from the observed mul- 
tiplicity, one should not expect a perfect quantitative 
agreement. In this sense the level of the agreement we 
obtained is quite acceptable and we would say that our 
results are consistent with the empirical values. 

A rough estimate of the average energy per gluon is 
(dE/drj) /{dN/drj) ~ 2.8 GeV. This value is greater than 
the largest energy anticipated in Ref. [ll], above which 
the CGC description might not make sense. This specu- 
lation is not quite correct, however. As is known and also 
clear in Fig. 21 the MV model is smoothly connected to 
the leading-order perturbative calculation with increas- 
ing k±, so that the spectrum has a perturbative tail in 
the high momentum region and the gluon energy is not 
bounded by the saturation scale. In contrast, the nonlin- 
ear effect of the CGC will be the most important for the 
gluons with momentum of order of Qsj 3'iid the average 
energy per gluon scales with Qs, which is indeed the case 
in our formulae. 

It is also intriguing to mention on the "formation 
time", which can roughly read from Fig. [8] for the a = 25 
case. We see that the curve gets flattened around mr ~ 
0.2. Since m = 0.36Qs, the "formation time" denoted 
as td which controls the convergence of t£{t) (see dis- 
cussions in the second paper of Ref. [l2|) is estimated as 
td — 0.56/Qs — 0.1 fm or g'^jiAT — 0.8. This is pretty 
short as compared with the expected thermalization time 
< 1 fm. This quantitative estimate of the "formation 
time" might deserve further investigation in the context 
of the early thermalization problem. 

The future applications, in addition to the early 
thermalization problem, cover various aspects of phe- 
nomenology as follows. 1) The rapidity dependence with 
changing jiA by Qs{x) gives the whole multiplicity distri- 
bution. 2) Finite size p^^^^\x±/) in the transverse plane 
yields the initial eccentricity at finite impact parameter. 
3) The energy-momentum tensor specifies the initial con- 
dition connected to the kinematic or hydrodynamic de- 
scription. 4) The (augmented) UV time evolution of the 
background fields is useful to attack the rapidity depen- 
dent instability on top of them [2^. What is addressed 
here is the starting point for all of these issues. We note 
that the point 2) is especially important because the ini- 
tial eccentricity obtained in the KLN model [lOl. [3^ is 
prevailing as the CGC estimate at present [40|, though 
the confirmation from the MV model is definitely neces- 
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sary [4]| . 

Finally, we close our discussions with a conclusion that 
our analytical estimates turn out to be consistent with 
the numerical simulations in spite of the difference in 
their model definitions [2^. This is an important check 
for the MV model foundation. 
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